Code for “Investigating Current Clinical Opinions in Stereo- Electroencephalography-informed Epilepsy Surgery”

This document provides the code that generated the results for the manuscript entitled “Investigating Current Clinical Opinions in Stereo- Electroencephalography-informed Epilepsy Surgery”. All code chunks and printed messages from computations are folded by default. One can click to unfold them to see the details. All code used for data wrangling, visualization, model fitting, model diagnostics, and model inference are included. A seed is provided for each model to maximize the numerical reproducibility of the MCMC model fitting. The original computation environment, including the operating system and software versions, is listed at the end of the document. The data are saved as SurveyData.xlsx after removing the country and gender information, one can load it as df (the resulting data frame of section Data wrangling) to run the rest of computations.

Library and utility function

PieDonut_new <- function(data, mapping, start = getOption("PieDonut.start",
    0), addPieLabel = TRUE, addDonutLabel = TRUE, showRatioDonut = TRUE,
    showRatioPie = TRUE, ratioByGroup = TRUE, showRatioThreshold = getOption("PieDonut.showRatioThreshold",
        0.02), labelposition = getOption("PieDonut.labelposition",
        2), labelpositionThreshold = 0.1, r0 = getOption("PieDonut.r0",
        0.3), r1 = getOption("PieDonut.r1", 1), r2 = getOption("PieDonut.r2",
        1.2), explode = NULL, selected = NULL, explodePos = 0.1,
    color = "white", pieAlpha = 0.8, donutAlpha = 1, maxx = NULL,
    showPieName = TRUE, showDonutName = FALSE, title = NULL,
    pieLabelSize = 4, donutLabelSize = 3, titlesize = 5, explodePie = TRUE,
    explodeDonut = FALSE, use.label = TRUE, use.labels = TRUE,
    family = getOption("PieDonut.family", ""), accuracy = 0.1) {
    (cols = colnames(data))
    if (use.labels)
        data = addLabelDf(data, mapping)
    count <- NULL
    if ("count" %in% names(mapping))
        count <- getMapping(mapping, "count")
    count
    pies <- donuts <- NULL
    (pies = getMapping(mapping, "pies"))
    if (is.null(pies))
        (pies = getMapping(mapping, "pie"))
    if (is.null(pies))
        (pies = getMapping(mapping, "x"))
    (donuts = getMapping(mapping, "donuts"))
    if (is.null(donuts))
        (donuts = getMapping(mapping, "donut"))
    if (is.null(donuts))
        (donuts = getMapping(mapping, "y"))
    if (!is.null(count)) {
        df <- data %>%
            group_by(.data[[pies]]) %>%
            dplyr::summarize(Freq = sum(.data[[count]]))
        df
    } else {
        df = data.frame(table(data[[pies]]))
    }
    colnames(df)[1] = pies
    df$end = cumsum(df$Freq)
    df$start = dplyr::lag(df$end)
    df$start[1] = 0
    total = sum(df$Freq)
    df$start1 = df$start * 2 * pi/total
    df$end1 = df$end * 2 * pi/total
    df$start1 = df$start1 + start
    df$end1 = df$end1 + start
    df$focus = 0
    if (explodePie)
        df$focus[explode] = explodePos
    df$mid = (df$start1 + df$end1)/2
    df$x = ifelse(df$focus == 0, 0, df$focus * sin(df$mid))
    df$y = ifelse(df$focus == 0, 0, df$focus * cos(df$mid))
    df$label = df[[pies]]
    df$ratio = df$Freq/sum(df$Freq)
    if (showRatioPie) {
        df$label = ifelse(df$ratio >= showRatioThreshold, paste0(df$label,
            "\n(", scales::percent(df$ratio, accuracy = accuracy),
            ")"), as.character(df$label))
    }
    df$labelx = (r0 + r1)/2 * sin(df$mid) + df$x
    df$labely = (r0 + r1)/2 * cos(df$mid) + df$y
    if (!is.factor(df[[pies]]))
        df[[pies]] <- factor(df[[pies]])
    df
    mainCol = gg_color_hue(nrow(df))
    df$radius = r1
    df$radius[df$focus != 0] = df$radius[df$focus != 0] + df$focus[df$focus !=
        0]
    df$hjust = ifelse((df$mid%%(2 * pi)) > pi, 1, 0)
    df$vjust = ifelse(((df$mid%%(2 * pi)) < (pi/2)) | (df$mid%%(2 *
        pi) > (pi * 3/2)), 0, 1)
    df$segx = df$radius * sin(df$mid)
    df$segy = df$radius * cos(df$mid)
    df$segxend = (df$radius + 0.05) * sin(df$mid)
    df$segyend = (df$radius + 0.05) * cos(df$mid)
    df
    if (!is.null(donuts)) {
        subColor = makeSubColor(mainCol, no = length(unique(data[[donuts]])))
        subColor
        data
        if (!is.null(count)) {
            df3 <- as.data.frame(data[c(donuts, pies, count)])
            colnames(df3) = c("donut", "pie", "Freq")
            df3
            df3 <- eval(parse(text = "complete(df3,donut,pie)"))
            df3$Freq[is.na(df3$Freq)] = 0
            if (!is.factor(df3[[1]]))
                df3[[1]] = factor(df3[[1]])
            if (!is.factor(df3[[2]]))
                df3[[2]] = factor(df3[[2]])
            df3 <- df3 %>%
                arrange(.data$pie, .data$donut)
            a <- df3 %>%
                spread(.data$pie, value = .data$Freq)
            a = as.data.frame(a)
            a
            rownames(a) = a[[1]]
            a = a[-1]
            a
            colnames(df3)[1:2] = c(donuts, pies)
        } else {
            df3 = data.frame(table(data[[donuts]], data[[pies]]),
                stringsAsFactors = FALSE)
            colnames(df3)[1:2] = c(donuts, pies)
            a = table(data[[donuts]], data[[pies]])
            a
        }
        a
        df3
        df3$group = rep(colSums(a), each = nrow(a))
        df3$pie = rep(1:ncol(a), each = nrow(a))
        total = sum(df3$Freq)
        total
        df3$ratio1 = df3$Freq/total
        df3
        if (ratioByGroup) {
            df3$ratio = scales::percent(df3$Freq/df3$group, accuracy = accuracy)
        } else {
            df3$ratio <- scales::percent(df3$ratio1, accuracy = accuracy)
        }
        df3$end = cumsum(df3$Freq)
        df3
        df3$start = dplyr::lag(df3$end)
        df3$start[1] = 0
        df3$start1 = df3$start * 2 * pi/total
        df3$end1 = df3$end * 2 * pi/total
        df3$start1 = df3$start1 + start
        df3$end1 = df3$end1 + start
        df3$mid = (df3$start1 + df3$end1)/2
        df3$focus = 0
        if (!is.null(selected)) {
            df3$focus[selected] = explodePos
        } else if (!is.null(explode)) {
            selected = c()
            for (i in 1:length(explode)) {
                start = 1 + nrow(a) * (explode[i] - 1)
                selected = c(selected, start:(start + nrow(a) -
                  1))
            }
            selected
            df3$focus[selected] = explodePos
        }
        df3
        df3$x = 0
        df3$y = 0
        df
        if (!is.null(explode)) {
            explode
            for (i in 1:length(explode)) {
                xpos = df$focus[explode[i]] * sin(df$mid[explode[i]])
                ypos = df$focus[explode[i]] * cos(df$mid[explode[i]])
                df3$x[df3$pie == explode[i]] = xpos
                df3$y[df3$pie == explode[i]] = ypos
            }
        }
        df3$no = 1:nrow(df3)
        df3$label = df3[[donuts]]
        if (showRatioDonut) {
            if (max(nchar(levels(df3$label))) <= 2)
                df3$label = paste0(df3$label, "(", df3$ratio,
                  ")") else df3$label = paste0(df3$label, "\n(", df3$ratio,
                ")")
        }
        df3$label[df3$ratio1 == 0] = ""
        df3$label[df3$ratio1 < showRatioThreshold] = ""
        df3$hjust = ifelse((df3$mid%%(2 * pi)) > pi, 1, 0)
        df3$vjust = ifelse(((df3$mid%%(2 * pi)) < (pi/2)) | (df3$mid%%(2 *
            pi) > (pi * 3/2)), 0, 1)
        df3$no = factor(df3$no)
        df3
        labelposition
        if (labelposition > 0) {
            df3$radius = r2
            if (explodeDonut)
                df3$radius[df3$focus != 0] = df3$radius[df3$focus !=
                  0] + df3$focus[df3$focus != 0]
            df3$segx = df3$radius * sin(df3$mid) + df3$x
            df3$segy = df3$radius * cos(df3$mid) + df3$y
            df3$segxend = (df3$radius + 0.05) * sin(df3$mid) +
                df3$x
            df3$segyend = (df3$radius + 0.05) * cos(df3$mid) +
                df3$y
            if (labelposition == 2)
                df3$radius = (r1 + r2)/2
            df3$labelx = (df3$radius) * sin(df3$mid) + df3$x
            df3$labely = (df3$radius) * cos(df3$mid) + df3$y
        } else {
            df3$radius = (r1 + r2)/2
            if (explodeDonut)
                df3$radius[df3$focus != 0] = df3$radius[df3$focus !=
                  0] + df3$focus[df3$focus != 0]
            df3$labelx = df3$radius * sin(df3$mid) + df3$x
            df3$labely = df3$radius * cos(df3$mid) + df3$y
        }
        df3$segx[df3$ratio1 == 0] = 0
        df3$segxend[df3$ratio1 == 0] = 0
        df3$segy[df3$ratio1 == 0] = 0
        df3$segyend[df3$ratio1 == 0] = 0
        if (labelposition == 0) {
            df3$segx[df3$ratio1 < showRatioThreshold] = 0
            df3$segxend[df3$ratio1 < showRatioThreshold] = 0
            df3$segy[df3$ratio1 < showRatioThreshold] = 0
            df3$segyend[df3$ratio1 < showRatioThreshold] = 0
        }
        df3
        del = which(df3$Freq == 0)
        del
        if (length(del) > 0)
            subColor <- subColor[-del]
        subColor
    }
    p <- ggplot() + theme_no_axes() + coord_fixed()
    if (is.null(maxx)) {
        r3 = r2 + 0.3
    } else {
        r3 = maxx
    }
    p1 <- p + geom_arc_bar(aes_string(x0 = "x", y0 = "y", r0 = as.character(r0),
        r = as.character(r1), start = "start1", end = "end1",
        fill = pies), alpha = pieAlpha, color = color, data = df) +
        transparent() + scale_fill_manual(values = mainCol) +
        xlim(r3 * c(-1, 1)) + ylim(r3 * c(-1, 1)) + guides(fill = FALSE)
    if ((labelposition == 1) & (is.null(donuts))) {
        p1 <- p1 + geom_segment(aes_string(x = "segx", y = "segy",
            xend = "segxend", yend = "segyend"), data = df) +
            geom_text(aes_string(x = "segxend", y = "segyend",
                label = "label", hjust = "hjust", vjust = "vjust"),
                size = pieLabelSize, data = df, family = family)
    } else if ((labelposition == 2) & (is.null(donuts))) {
        p1 <- p1 + geom_segment(aes_string(x = "segx", y = "segy",
            xend = "segxend", yend = "segyend"), data = df[df$ratio <
            labelpositionThreshold, ]) + geom_text(aes_string(x = "segxend",
            y = "segyend", label = "label", hjust = "hjust",
            vjust = "vjust"), size = pieLabelSize, data = df[df$ratio <
            labelpositionThreshold, ], family = family) + geom_text(aes_string(x = "labelx",
            y = "labely", label = "label"), size = pieLabelSize,
            data = df[df$ratio >= labelpositionThreshold, ],
            family = family)
    } else {
        p1 <- p1 + geom_text(aes_string(x = "labelx", y = "labely",
            label = "label"), size = pieLabelSize, data = df,
            family = family)
    }
    if (showPieName)
        p1 <- p1 + annotate("text", x = 0, y = 0, label = pies,
            size = titlesize, family = family)
    p1 <- p1 + theme(text = element_text(family = family))
    if (!is.null(donuts)) {
        if (explodeDonut) {
            p3 <- p + geom_arc_bar(aes_string(x0 = "x", y0 = "y",
                r0 = as.character(r1), r = as.character(r2),
                start = "start1", end = "end1", fill = "no",
                explode = "focus"), alpha = donutAlpha, color = color,
                data = df3)
        } else {
            p3 <- p + geom_arc_bar(aes_string(x0 = "x", y0 = "y",
                r0 = as.character(r1), r = as.character(r2),
                start = "start1", end = "end1", fill = "no"),
                alpha = donutAlpha, color = color, data = df3)
        }
        p3 <- p3 + transparent() + scale_fill_manual(values = subColor) +
            xlim(r3 * c(-1, 1)) + ylim(r3 * c(-1, 1)) + guides(fill = FALSE)
        p3
        if (labelposition == 1) {
            p3 <- p3 + geom_segment(aes_string(x = "segx", y = "segy",
                xend = "segxend", yend = "segyend"), data = df3) +
                geom_text(aes_string(x = "segxend", y = "segyend",
                  label = "label", hjust = "hjust", vjust = "vjust"),
                  size = donutLabelSize, data = df3, family = family)
        } else if (labelposition == 0) {
            p3 <- p3 + geom_text(aes_string(x = "labelx", y = "labely",
                label = "label"), size = donutLabelSize, data = df3,
                family = family)
        } else {
            p3 <- p3 + geom_segment(aes_string(x = "segx", y = "segy",
                xend = "segxend", yend = "segyend"), data = df3[df3$ratio1 <
                labelpositionThreshold, ]) + geom_text(aes_string(x = "segxend",
                y = "segyend", label = "label", hjust = "hjust",
                vjust = "vjust"), size = donutLabelSize, data = df3[df3$ratio1 <
                labelpositionThreshold, ], family = family) +
                geom_text(aes_string(x = "labelx", y = "labely",
                  label = "label"), size = donutLabelSize, data = df3[df3$ratio1 >=
                  labelpositionThreshold, ], family = family)
        }
        if (!is.null(title))
            p3 <- p3 + annotate("text", x = 0, y = r3, label = title,
                size = titlesize, family = family) else if (showDonutName)
            p3 <- p3 + annotate("text", x = (-1) * r3, y = r3,
                label = donuts, hjust = 0, size = titlesize,
                family = family)
        p3 <- p3 + theme(text = element_text(family = family))
        grid.newpage()
        print(p1, vp = viewport(height = 1, width = 1))
        print(p3, vp = viewport(height = 1, width = 1))
    } else {
        p1
    }
}


pie_plot <- function(df, met, title,
    legend.title) {
    df_pie <- df %>%
        group_by(!!sym(met)) %>%
        summarize(n = n()) %>%
        ungroup()

    tmp <- ggpiestats(df %>%
        filter(!!sym(met) != "Uncertain",
            !!sym(met) != "Unavailable",
            !!sym(met) != "Others", !!sym(met) !=
                "No choice"), !!sym(met))

    start <- case_when(met == "DCFS" ~
        pi * 3/4, met == "SOZname" ~
        pi * 3/4, met == "Continent" ~
        pi * 4/5, met == "LVFA" ~ pi *
        16/15, met == "Gender" ~ pi *
        4/5, met == "Age" ~ pi * 4/5,
        met == "SEEG" ~ pi * 4/5, .default = pi *
            3/4)
    if (met == "LVFA") {
        PieDonut_new(df_pie, aes(!!sym(met),
            count = n), r0 = 0.75, r1 = 1.3,
            r2 = 1.3, start = start,
            showRatioThreshold = 0.001,
            pieLabelSize = 10, donutLabelSize = 8,
            titlesize = 0, use.label = F,
            use.labels = F) + annotate("text",
            x = 0, y = 0, label = title,
            size = 13, lineheight = 0.8) +
            scale_color_manual(values = c("#374E55FF",
                "#DF8F44FF", "#00A1D5FF",
                "#B24745FF", "#6A6599FF",
                "#79AF97FF", "#80796BFF")) +
            scale_fill_manual(values = c("#374E55FF",
                "#DF8F44FF", "#00A1D5FF",
                "#B24745FF", "#6A6599FF",
                "#79AF97FF", "#80796BFF")) +
            theme(plot.title = element_blank(),
                plot.caption = element_text(size = 18)) ->
                p
    } else {
        PieDonut_new(df_pie, aes(!!sym(met),
            count = n), r0 = 0.75, r1 = 1.3,
            r2 = 1.3, start = start,
            showRatioThreshold = 0.001,
            pieLabelSize = 10, donutLabelSize = 8,
            titlesize = 0, use.label = F,
            use.labels = F) + annotate("text",
            x = 0, y = 0, label = title,
            size = 13, lineheight = 0.8) +
            scale_color_jama() + scale_fill_jama() +
            theme(plot.title = element_blank(),
                plot.caption = element_text(size = 18)) ->
                p
    }

    p_file <- tempfile(tmpdir = "./fig",
        fileext = ".png")
    agg_png(p_file, width = 8.5, height = 8.5,
        units = "in", res = 100)
    suppressWarnings(print(p))
    invisible(dev.off())

    return(list(p_file = p_file, p = p))
}


inference_plot_categorical <- function(draws, level_plt, title,
    height = 8) {
    drawsa <- draws %>%
        pivot_longer(cols = everything(), values_to = "value",
            names_to = "para") %>%
        mutate(para = factor(para, levels = level_plt))

    lims <- drawsa %>%
        group_by(para) %>%
        summarise(q1 = quantile(value, 0.001), q2 = quantile(value,
            0.999)) %>%
        ungroup() %>%
        summarise(limlow = min(q1), limhigh = max(q2))

    drawsa %>%
        ggplot(aes(x = value, y = para)) + stat_dotsinterval(quantiles = 100,
        point_interval = "mode_hdi", layout = "weave") + scale_x_continuous(breaks = scales::pretty_breaks(n = 5),
        limits = c(lims$limlow, lims$limhigh)) + geom_vline(xintercept = 0,
        linetype = "dashed") + labs(title = title, x = "Proportion") +
        theme(legend.position = "none", plot.margin = margin(5.5,
            50, 5.5, 5.5, "pt"), axis.title.y = element_blank(),
            axis.text.y = element_text(hjust = 0), panel.border = element_blank(),
            panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
            axis.line = element_blank(), panel.background = element_blank(),
            text = element_text(size = 38, face = "plain", color = "black"),
            axis.text = element_text(color = "black", size = 38),
            plot.title.position = "plot", plot.title = element_text(size = 32),
            plot.subtitle = element_text(size = 20)) -> p_tempa

    drawsb <- draws %>%
        mutate_at(vars(-Population), ~. - Population) %>%
        select(-Population) %>%
        pivot_longer(cols = everything(), values_to = "value",
            names_to = "para") %>%
        mutate(para = factor(para, levels = level_plt))

    lims <- drawsb %>%
        group_by(para) %>%
        summarise(q1 = quantile(value, 0.001), q2 = quantile(value,
            0.999)) %>%
        ungroup() %>%
        summarise(limlow = min(q1), limhigh = max(q2))

    pd <- draws %>%
        mutate_at(vars(-Population), ~. - Population) %>%
        select(-Population) %>%
        mutate_all(~. > 0) %>%
        summarize_all(~mean(.))

    drawsb %>%
        ggplot(aes(x = value, y = para, fill = after_stat(x >
            0))) + stat_dotsinterval(quantiles = 100, point_interval = "mode_hdi",
        layout = "weave") + scale_fill_manual(values = c("#DF8F44FF",
        "#00A1D5FF")) + scale_x_continuous(breaks = scales::pretty_breaks(n = 5),
        limits = c(lims$limlow, lims$limhigh)) + geom_vline(xintercept = 0,
        linetype = "dashed") + annotate("text", x = 0, y = colnames(pd),
        label = scales::percent(unlist(pd), accuracy = 0.1),
        size = 10, vjust = -2.5, hjust = -1) + annotate("text",
        x = 0, y = colnames(pd), label = scales::percent(1 -
            unlist(pd), accuracy = 0.1), size = 10, vjust = -2.5,
        hjust = 2) + labs(title = " Deviation from population",
        x = "Proportion") + theme(legend.position = "none", plot.margin = margin(5.5,
        50, 5.5, 5.5, "pt"), axis.title.y = element_blank(),
        axis.text.y = element_text(hjust = 0), panel.border = element_blank(),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.line = element_blank(), panel.background = element_blank(),
        text = element_text(size = 38, face = "plain", color = "black"),
        axis.text = element_text(color = "black", size = 38),
        plot.title.position = "plot", plot.title = element_text(size = 32),
        plot.subtitle = element_text(size = 20)) -> p_tempb

    p_temp <- ggarrange(p_tempa, p_tempb, nrow = 1, ncol = 2)

    p_temp_file <- tempfile(tmpdir = "./fig", fileext = ".png")
    agg_png(p_temp_file, width = 20, height = height, units = "in",
        res = 100)
    suppressWarnings(print(p_temp))
    invisible(dev.off())

    return(list(p_temp = p_temp, p_temp_file = p_temp_file))
}

inference_plot_categorical_one <- function(draws, level_plt,
    title, height = 8) {
    color_code <- propotion %>%
        mutate(colorcode = case_when(Continent == "America" ~
            "#DF8F44FF", Continent == "Asia Pacific" ~ "#00A1D5FF",
            Continent == "Europe" ~ "#374E55FF", TRUE ~ as.character(Continent))) %>%
        rename(para = Country) %>%
        select(para, colorcode)

    drawsa <- draws %>%
        pivot_longer(cols = everything(), values_to = "value",
            names_to = "para") %>%
        left_join(color_code, by = c("para")) %>%
        mutate(para = factor(para, levels = level_plt))

    lims <- drawsa %>%
        group_by(para) %>%
        summarise(q1 = quantile(value, 0.001), q2 = quantile(value,
            0.999)) %>%
        ungroup() %>%
        summarise(limlow = min(q1), limhigh = max(q2))

    drawsa %>%
        ggplot(aes(x = value, y = para)) + stat_dotsinterval(quantiles = 100,
        point_interval = "mode_hdi", layout = "weave") + scale_x_continuous(breaks = scales::pretty_breaks(n = 5),
        limits = c(lims$limlow, lims$limhigh)) + geom_vline(xintercept = as.numeric(mode_hdi(draws$Population)[1,
        1]), linetype = "dashed") + labs(title = title, x = "Proportion") +
        theme(legend.position = "none", plot.margin = margin(5.5,
            50, 5.5, 5.5, "pt"), axis.title.y = element_blank(),
            axis.text.y = element_text(hjust = 0), panel.border = element_blank(),
            panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
            axis.line = element_blank(), panel.background = element_blank(),
            text = element_text(size = 38, face = "plain", color = "black"),
            axis.text = element_text(color = "black", size = 38),
            plot.title.position = "plot", plot.title = element_text(size = 32),
            plot.subtitle = element_text(size = 20)) -> p_tempa

    p_temp_file <- tempfile(tmpdir = "./fig", fileext = ".png")
    agg_png(p_temp_file, width = 10, height = height, units = "in",
        res = 100)
    suppressWarnings(print(p_tempa))
    invisible(dev.off())

    return(list(p_temp = p_tempa, p_temp_file = p_temp_file))
}

plot_hist <- function(df,
    ref, title) {

    percentile_intervals <- hdi(df)

    if (mean(df > ref) >
        0.5) {
        hjust_tmp <- -0.01
        color_tmp <- "#00A1D5FF"
        probtext <- "\nProb(MRP > Observed) = "
        prop_tmp <- mean(df >
            ref)
    } else {
        hjust_tmp <- 1.01
        color_tmp <- "#DF8F44FF"
        probtext <- "\nProb(MRP < Observed) = "
        prop_tmp <- mean(df <
            ref)
    }

    data.frame(x = df) %>%
        ggplot(aes(x = x,
            fill = after_stat(x >
                ref))) +
        stat_dotsinterval(quantiles = 100,
            point_interval = "mode_hdi",
            layout = "weave",
            interval_size_range = c(1.2,
                2.8)) +
        scale_fill_manual(values = c("#DF8F44FF",
            "#00A1D5FF")) +
        labs(title = title,
            x = "Proportion") +
        geom_vline(xintercept = ref,
            col = color_tmp,
            size = 1) +
        annotate("text",
            x = round(as.numeric(mode_hdi(df)[1]),
                digits = 3),
            y = 0.1, label = paste0("MRP: ",
                sprintf("%0.3f",
                  round(as.numeric(mode_hdi(df)[1]),
                    digits = 3))),
            size = 12, hjust = 0.5,
            col = "black") +
        annotate("text",
            x = ref, y = 1.05,
            label = paste0("Observed proportion: ",
                sprintf("%0.3f",
                  round(ref,
                    digits = 3)),
                probtext,
                scales::percent(prop_tmp,
                  0.1)),
            size = 10, hjust = hjust_tmp,
            color = color_tmp) +
        scale_x_continuous(breaks = breaks_pretty(4),
            expand = c(0.01,
                0.01)) +
        scale_y_continuous(limits = c(NA,
            1.1)) + labs(subtitle = " ") +
        theme_classic() +
        theme(legend.position = "none",
            plot.margin = margin(5.5,
                20, 5.5,
                20, "pt"),
            axis.title.y = element_blank(),
            axis.ticks.y = element_blank(),
            axis.text.y = element_blank(),
            panel.border = element_blank(),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            axis.line = element_blank(),
            panel.background = element_blank(),
            text = element_text(size = 38,
                face = "plain",
                color = "black"),
            axis.text = element_text(color = "black",
                size = 38),
            plot.title.position = "plot",
            plot.title = element_text(size = 36),
            plot.subtitle = element_text(size = 28)) ->
            p_temp

    p_temp_file <- tempfile(tmpdir = "./fig",
        fileext = ".png")
    agg_png(p_temp_file,
        width = 12, height = 7.1,
        units = "in", res = 100)
    suppressWarnings(print(p_temp))
    invisible(dev.off())

    return(list(p_temp = p_temp,
        p_temp_file = p_temp_file))
}

Data wrangling

Load data and mutate variable values and levels.

# load data
propotion <- read_excel("./Survey.xlsx", sheet = "Sheet1") %>% 
  rename(Continent = `Continent ID`,
         Responses = `No of responses`,
         Centers = `No of SEEG centers`,
         Country = Country) %>% 
  select(Country, Continent, Responses, Centers) %>% 
  mutate(Country = str_replace_all(str_to_title(Country), "'", "")) %>% 
  mutate(Continent = case_when(
    Continent %in% c(1, 5) ~ "America",
    Continent == 2 ~ "Europe",
    Continent %in% c(3, 4) ~ "Asia Pacific",
    is.na(Continent) ~ "Unavailable",
    TRUE ~ as.character(Continent)
  )) %>% 
  filter(!Country %in% c("Latvia", "Mexico", "Poland", "South Korea"))

df <- read_excel("./Survey.xlsx", sheet = "Sheet2") %>% 
  rename(ResponseTime = `Response time`,
         Age = `Age group`,
         SEEG = `SEEEG number at their center`,
         Surgery = `% SEEG resulting in surgery`,
         DCFS = `DC/Fs`,
         LVFA = `LVFA def`,
         EZ = `EZ def`,
         SOZ = `SOZ def`,
         SZEZ = `SZ/EZ confusion`,
         SOZname = `SOZ name`) %>% 
  mutate(Age = case_when(
    Age == 0 ~ "Unavailable",
    Age == 1 ~ "25 - 35",
    Age == 2 ~ "35 - 50",
    Age == 3 ~ "> 50",
    TRUE ~ as.character(Age)
  )) %>% 
  mutate(Gender = case_when(
    Gender == 0 ~ "Unavailable",
    Gender == 1 ~ "Male",
    Gender == 2 ~ "Female",
    TRUE ~ as.character(Gender)
  )) %>%
  mutate(Experience = case_when(
    Experience == 0 ~ "Unavailable",
    Experience == 1 ~ "Low",
    Experience == 2 ~ "Moderate",
    Experience == 3 ~ "Moderate",
    Experience == 4 ~ "High",
    Experience == 5 ~ "High",
    TRUE ~ as.character(Experience)
  )) %>% 
  mutate(SEEG = case_when(
    SEEG == 0 ~ "Unavailable",
    SEEG == 1 ~ "< 10",
    SEEG == 2 ~ "10 - 25",
    SEEG == 3 ~ "25 - 40",
    SEEG == 4 ~ "> 40",
    TRUE ~ as.character(SEEG)
  )) %>% 
  mutate(Surgery = case_when(
    Surgery == 0 ~ "Unavailable",
    Surgery == 1 ~ "> 80%",
    Surgery == 2 ~ "50% - 80%",
    Surgery == 3 ~ "< 50%",
    TRUE ~ as.character(Surgery)
  )) %>% 
  mutate(DCFS = case_when(
    DCFS == 0 ~ "Uncertain",
    DCFS == 1 ~ "DC | Fs<1kHz",
    DCFS == 2 ~ "noDC | Fs>1kHz",
    DCFS == 3 ~ "Both",
    DCFS == 4 ~ "Neither",
    TRUE ~ as.character(DCFS)
  )) %>% 
  mutate(LVFA = case_when(
    LVFA == 0 ~ "Unavailable",
    LVFA == 1 ~ "13Hz",
    LVFA == 2 ~ "30Hz",
    LVFA == 3 ~ "60Hz",
    LVFA == 4 ~ "100Hz",
    LVFA == 5 ~ "Not use",
    LVFA == 6 ~ "Others",
    TRUE ~ as.character(LVFA)
  )) %>% 
  mutate(EZ = case_when(
    EZ == 0 ~ "Others",
    EZ == 1 ~ "Lüders",
    EZ == 2 ~ "Bancaud",
    TRUE ~ as.character(EZ)
  )) %>% 
  mutate(SOZ = case_when(
    SOZ == 0 ~ "Others",
    SOZ == 1 ~ "Onset",
    SOZ == 2 ~ "Onset + P1",
    TRUE ~ as.character(SOZ)
  )) %>% 
  mutate(SZEZ = case_when(
    SZEZ == 0 ~ "No",
    SZEZ == 1 ~ "Yes",
    SZEZ == 2 ~ "Maybe",
    TRUE ~ as.character(SZEZ)
  )) %>% 
  mutate(SOZname = case_when(
    SOZname == 0 ~ "Others",
    SOZname == 1 ~ "HEZ",
    SOZname == 2 ~ "PEZ",
    SOZname == 3 ~ "EZR",
    TRUE ~ as.character(SOZname)
  ))

df <- df %>% 
  mutate(Age = factor(Age, levels = c("Unavailable", "25 - 35", "35 - 50", "> 50")),
         Gender = factor(Gender, levels = c("Unavailable", "Male", "Female")),
         DCFS = factor(DCFS, levels = 
                         c("Neither", "DC | Fs<1kHz", "noDC | Fs>1kHz", "Both", "Uncertain")),
         Experience = factor(Experience, levels = c("Unavailable", "Low", "Moderate", "High")),
         SEEG = factor(SEEG, levels = c("Unavailable", "< 10", "10 - 25", "25 - 40", "> 40")),
         Surgery = factor(Surgery, levels = c("Unavailable", "< 50%", "50% - 80%", "> 80%")),
         EZ = factor(EZ, levels = c("Others", "Lüders", "Bancaud")),
         SOZ = factor(SOZ, levels = c("Others", "Onset + P1", "Onset")),
         SZEZ = factor(SZEZ, levels = c("Maybe", "No", "Yes")),
         SOZname = factor(SOZname, levels = c("Others", "HEZ", "PEZ", "EZR")),
         LVFA = factor(LVFA, levels = c("Unavailable", "13Hz", "30Hz",  
                                       "60Hz", "100Hz", "Others", "Not use")))

df <- df %>% 
  mutate(Country = str_replace_all(str_to_title(Country), "'", "")) %>% 
  mutate(Country = ifelse(Country == "", "Unavailable", Country)) %>% 
  left_join(propotion, by = c("Country")) %>% 
  mutate(Continent = ifelse(is.na(Continent), "Unavailable", Continent)) %>% 
  mutate(Country = factor(Country),
         Continent = factor(Continent, levels = c("Unavailable", "America", "Asia Pacific", "Europe")))

feature <- read_excel("./Survey.xlsx", sheet = "feature")

feature_names <- tibble(name = colnames(feature),
                        index = paste0("feature", seq(1, 14)))

feature <- feature %>% 
  rename_all(~feature_names$index) %>% 
  mutate_all(~ifelse(. == 0, "No", "Yes")) %>% 
  mutate_all(~factor(., levels = c("No", "Yes")))

df <- df %>% 
  add_column(feature)

Distributions of the responses

Figure 2

ptmp_file <- pie_plot(df, "Gender", title = "Gender",
                      legend.title = "Gender")
knitr::include_graphics(ptmp_file$p_file)

ptmp_file <- pie_plot(df, "Age", title = "Age\n(year)",
                      legend.title = "Age")
knitr::include_graphics(ptmp_file$p_file)

ptmp_file <- pie_plot(df, "DCFS", title = "DC/Fs",
                      legend.title = "DCFS")
knitr::include_graphics(ptmp_file$p_file)

ptmp_file <- pie_plot(df, "Experience", title = "Experience",
                      legend.title = "Experience")
knitr::include_graphics(ptmp_file$p_file)

ptmp_file <- pie_plot(df, "SEEG", title = "SEEG\n(volume/year)",
                      legend.title = "SEEG")
knitr::include_graphics(ptmp_file$p_file)

ptmp_file <- pie_plot(df, "Surgery", title = "Post-SEEG\n(surgery/year)",
                      legend.title = "Surgery")
knitr::include_graphics(ptmp_file$p_file)

ptmp_file <- pie_plot(df, "EZ", title = "EZ",
                      legend.title = "EZ")
knitr::include_graphics(ptmp_file$p_file)

ptmp_file <- pie_plot(df, "SOZ", title = "SOZ",
                      legend.title = "SOZ")
knitr::include_graphics(ptmp_file$p_file)

ptmp_file <- pie_plot(df, "Continent", title = "Geographic\nArea",
                      legend.title = "Continent")
knitr::include_graphics(ptmp_file$p_file)

ptmp_file <- pie_plot(df %>% 
                         mutate(LVFA = as.character(LVFA)) %>% 
                         mutate(LVFA = ifelse(LVFA == "Not use", "Do not use", LVFA)) %>% 
                         mutate(LVFA = factor(LVFA, levels = c("Unavailable", "13Hz", 
                                                                "30Hz", "60Hz", "Do not use", "100Hz", 
                                                                "Others"))), "LVFA", title = "LVFA",
                      legend.title = "LVFA")
knitr::include_graphics(ptmp_file$p_file)

ptmp_file <- pie_plot(df, "SZEZ", title = "SOZ/EZ\nconfusion",
                      legend.title = "SZEZ")
knitr::include_graphics(ptmp_file$p_file)

ptmp_file <- pie_plot(df, "SOZname", title = "SOZ\nterminology",
                      legend.title = "SOZname")
knitr::include_graphics(ptmp_file$p_file)

EZ definition across experience

Figure 3B

df_complete <- df %>% 
  filter(EZ != "Others",
         Experience != "Unavailable") %>% 
  mutate_if(is.factor, droplevels)

mod_fit <- brm(
  data = df_complete, family = bernoulli(link = "logit"),
  EZ ~ 1 + (1 | Experience),
  prior = c(prior(normal(0, 1.5), class = Intercept), 
            prior(normal(0, 1), class = sd)), 
  control = list(adapt_delta = 0.99),
  iter = 3500, warmup = 1000, chains = 4, cores = 4, seed = 20230918)

mod_fit

Model fitting results and diagnostics

 Family: bernoulli 
  Links: mu = logit 
Formula: EZ ~ 1 + (1 | Experience) 
   Data: df_complete (Number of observations: 304) 
  Draws: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
         total post-warmup draws = 10000

Group-Level Effects: 
~Experience (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.56      0.39     0.06     1.55 1.00     1852     2271

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    -0.71      0.40    -1.43     0.24 1.00     2004     1974

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

fdraws <- as_draws_df(mod_fit) %>% 
  as_tibble() %>% 
  select(-starts_with(".")) %>% 
  select(-starts_with("lp")) %>% 
  summarise(Population = inv_logit_scaled(b_Intercept),
            `Low` = inv_logit_scaled(b_Intercept + `r_Experience[Low,Intercept]`),
            `Moderate` = inv_logit_scaled(b_Intercept + `r_Experience[Moderate,Intercept]`),
            `High` = inv_logit_scaled(b_Intercept + `r_Experience[High,Intercept]`)) 


level_plt <- c("Population", "Low", "Moderate", "High")

p_temp <- inference_plot_categorical(fdraws, level_plt, paste0(tail(levels(df_complete$EZ), 1), " by Experience"))

knitr::include_graphics(p_temp$p_temp_file)

EZ definition across geographic area

Figure 3C and D

df_complete <- df %>% 
  filter(EZ != "Others",
         Country != "Unavailable") %>% 
  mutate_if(is.factor, droplevels)

mod_fit <- brm(
  data = df_complete, family = bernoulli(link = "logit"),
  EZ ~ 1 + (1 | Continent) + (1 | Country),
  prior = c(prior(normal(0, 1.5), class = Intercept), 
            prior(normal(0, 1), class = sd)), 
  control = list(adapt_delta = 0.99),
  iter = 3500, warmup = 1000, chains = 4, cores = 4, seed = 20230918)

mod_fit

Model fitting results and diagnostics

 Family: bernoulli 
  Links: mu = logit 
Formula: EZ ~ 1 + (1 | Continent) + (1 | Country) 
   Data: df_complete (Number of observations: 312) 
  Draws: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
         total post-warmup draws = 10000

Group-Level Effects: 
~Continent (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.54      0.46     0.02     1.72 1.00     3409     4520

~Country (Number of levels: 35) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     1.51      0.36     0.87     2.30 1.00     4591     6243

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    -1.27      0.54    -2.31    -0.15 1.00     4686     5094

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

fdraws <- as_draws_df(mod_fit) %>% 
  as_tibble() %>% 
  select(-starts_with(".")) %>% 
  select(-starts_with("lp")) %>% 
  mutate(Population = inv_logit_scaled(b_Intercept),
         "Asia Pacific" = inv_logit_scaled(b_Intercept + `r_Continent[Asia.Pacific,Intercept]`),
         "Europe" = inv_logit_scaled(b_Intercept + `r_Continent[Europe,Intercept]`),
         "America" = inv_logit_scaled(b_Intercept + `r_Continent[America,Intercept]`)) 

for (country in propotion$Country) {
  if (propotion$Responses[propotion$Country == country] >= 4){
      fdraws[[country]] <- inv_logit_scaled(fdraws[["b_Intercept"]] + 
                                          fdraws[[paste0("r_Continent[", 
                                                         str_replace(propotion$Continent[propotion$Country == country], " ", "."), 
                                                         ",Intercept]")]] + 
                                          fdraws[[paste0("r_Country[", 
                                                         str_replace(country, " ", "."), 
                                                         ",Intercept]")]])
  }

}

fdraws <- fdraws %>% 
  select(-starts_with(c("b_", "r_", "sd_")))

level_plt <- c("Population", "Asia Pacific", "Europe", "America")

p_temp <- inference_plot_categorical(fdraws %>% select("Population", "Asia Pacific", 
                                                       "Europe", "America"), 
                                     level_plt, paste0(tail(levels(df_complete$EZ), 1), " by Geographic area"))

knitr::include_graphics(p_temp$p_temp_file)

level_plt <- propotion %>% 
  arrange(desc(Continent), desc(Country)) %>% 
  filter(Country %in% colnames(fdraws %>% select(-c("Asia Pacific", "Europe", "America")))) %>% 
  pull(Country)


level_plt <- c("Population", level_plt)

p_temp <- inference_plot_categorical_one(fdraws %>% select(-c("Asia Pacific", "Europe", "America")), level_plt,
                                     " Country of employment",
                                     height = 30)
knitr::include_graphics(p_temp$p_temp_file)

SOZ defintion across experience

Figure S14

df_complete <- df %>% 
  filter(SOZ != "Others",
         Experience != "Unavailable") %>% 
  mutate_if(is.factor, droplevels)

mod_fit <- brm(
  data = df_complete, family = bernoulli(link = "logit"),
  SOZ ~ 1 + (1 | Experience),
  prior = c(prior(normal(0, 1.5), class = Intercept), 
            prior(normal(0, 1), class = sd)), 
  control = list(adapt_delta = 0.99),
  iter = 3500, warmup = 1000, chains = 4, cores = 4, seed = 20230918)

mod_fit

Model fitting results and diagnostics

 Family: bernoulli 
  Links: mu = logit 
Formula: SOZ ~ 1 + (1 | Experience) 
   Data: df_complete (Number of observations: 301) 
  Draws: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
         total post-warmup draws = 10000

Group-Level Effects: 
~Experience (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.29      0.29     0.01     1.07 1.00     2183     3147

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     0.30      0.25    -0.20     0.83 1.00     2578     2464

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

fdraws <- as_draws_df(mod_fit) %>% 
  as_tibble() %>% 
  select(-starts_with(".")) %>% 
  select(-starts_with("lp")) %>% 
  summarise(Population = inv_logit_scaled(b_Intercept),
            `Low` = inv_logit_scaled(b_Intercept + `r_Experience[Low,Intercept]`),
            `Moderate` = inv_logit_scaled(b_Intercept + `r_Experience[Moderate,Intercept]`),
            `High` = inv_logit_scaled(b_Intercept + `r_Experience[High,Intercept]`)) 

level_plt <- c("Population", "Low", "Moderate", "High")

p_temp <- inference_plot_categorical(fdraws, level_plt, paste0(tail(levels(df_complete$SOZ), 1), " by Experience"))

knitr::include_graphics(p_temp$p_temp_file)

SOZ defintion across geographic area

Figure S14

df_complete <- df %>% 
  filter(SOZ != "Others",
         Country != "Unavailable") %>% 
  mutate_if(is.factor, droplevels)

mod_fit <- brm(
  data = df_complete, family = bernoulli(link = "logit"),
  SOZ ~ 1 + (1 | Continent) + (1 | Country),
  prior = c(prior(normal(0, 1.5), class = Intercept), 
            prior(normal(0, 1), class = sd)), 
  control = list(adapt_delta = 0.99),
  iter = 3500, warmup = 1000, chains = 4, cores = 4, seed = 20230918)

mod_fit

Model fitting results and diagnostics

 Family: bernoulli 
  Links: mu = logit 
Formula: SOZ ~ 1 + (1 | Continent) + (1 | Country) 
   Data: df_complete (Number of observations: 309) 
  Draws: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
         total post-warmup draws = 10000

Group-Level Effects: 
~Continent (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.35      0.34     0.01     1.30 1.00     3971     4188

~Country (Number of levels: 35) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.63      0.26     0.12     1.17 1.00     2647     2129

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     0.33      0.33    -0.36     0.98 1.00     4319     4333

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

fdraws <- as_draws_df(mod_fit) %>% 
  as_tibble() %>% 
  select(-starts_with(".")) %>% 
  select(-starts_with("lp")) %>% 
  mutate(Population = inv_logit_scaled(b_Intercept),
         "Asia Pacific" = inv_logit_scaled(b_Intercept + `r_Continent[Asia.Pacific,Intercept]`),
         "Europe" = inv_logit_scaled(b_Intercept + `r_Continent[Europe,Intercept]`),
         "America" = inv_logit_scaled(b_Intercept + `r_Continent[America,Intercept]`)) 

for (country in propotion$Country) {
  if (propotion$Responses[propotion$Country == country] >= 4){
    fdraws[[country]] <- inv_logit_scaled(fdraws[["b_Intercept"]] + 
                                            fdraws[[paste0("r_Continent[", 
                                                           str_replace(propotion$Continent[propotion$Country == country], " ", "."), 
                                                           ",Intercept]")]] + 
                                            fdraws[[paste0("r_Country[", 
                                                           str_replace(country, " ", "."), 
                                                           ",Intercept]")]])
  }
  
}

fdraws <- fdraws %>% 
  select(-starts_with(c("b_", "r_", "sd_")))

level_plt <- c("Population", "Asia Pacific", "Europe", "America")

p_temp <- inference_plot_categorical(fdraws %>% select("Population", "Asia Pacific", "Europe", "America"), level_plt, paste0(tail(levels(df_complete$SOZ), 1), " by Geographic area"))

knitr::include_graphics(p_temp$p_temp_file)

Multivariate and multilevel regression and post-stratification (MRP) analysis

EZ defintion

df_complete <- df %>% 
  filter(EZ != "Others",
         Country != "Unavailable",
         Experience != "Unavailable") %>% 
  mutate_if(is.factor, droplevels)

mod_fit <- brm(
  data = df_complete, family = bernoulli(link = "logit"),
  EZ ~ 1 + (1 | Continent) + (1 | Country) + (1 | Experience),
  prior = c(prior(normal(0, 1.5), class = Intercept), 
            prior(normal(0, 1), class = sd)), 
  control = list(adapt_delta = 0.99),
  iter = 3500, warmup = 1000, chains = 4, cores = 4, seed = 20230918)

mod_fit

Model fitting results and diagnostics

 Family: bernoulli 
  Links: mu = logit 
Formula: EZ ~ 1 + (1 | Continent) + (1 | Country) + (1 | Experience) 
   Data: df_complete (Number of observations: 299) 
  Draws: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
         total post-warmup draws = 10000

Group-Level Effects: 
~Continent (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.55      0.45     0.02     1.69 1.00     4033     4859

~Country (Number of levels: 34) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     1.42      0.37     0.78     2.22 1.00     4393     5832

~Experience (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.50      0.39     0.03     1.56 1.00     4215     3621

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    -1.15      0.60    -2.31     0.10 1.00     5976     6234

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

fdraws <- as_draws_df(mod_fit) %>% 
  as_tibble() %>% 
  select(-starts_with(".")) %>% 
  select(-starts_with("lp")) %>% 
  summarise(Population = inv_logit_scaled(b_Intercept),
            `Low` = inv_logit_scaled(b_Intercept + `r_Experience[Low,Intercept]`),
            `Moderate` = inv_logit_scaled(b_Intercept + `r_Experience[Moderate,Intercept]`),
            `High` = inv_logit_scaled(b_Intercept + `r_Experience[High,Intercept]`)) 

level_plt <- c("Population", "Low", "Moderate", "High")

p_temp <- inference_plot_categorical(fdraws, level_plt, paste0(tail(levels(df_complete$EZ), 1), " by Experience"))

knitr::include_graphics(p_temp$p_temp_file)

fdraws <- as_draws_df(mod_fit) %>% 
  as_tibble() %>% 
  select(-starts_with(".")) %>% 
  select(-starts_with("lp")) %>% 
  mutate(Population = inv_logit_scaled(b_Intercept),
         "Asia Pacific" = inv_logit_scaled(b_Intercept + `r_Continent[Asia.Pacific,Intercept]`),
         "Europe" = inv_logit_scaled(b_Intercept + `r_Continent[Europe,Intercept]`),
         "America" = inv_logit_scaled(b_Intercept + `r_Continent[America,Intercept]`)) 

fdraws <- fdraws %>% 
  select(-starts_with(c("b_", "r_", "sd_")))

level_plt <- c("Population", "Asia Pacific", "Europe", "America")

p_temp <- inference_plot_categorical(fdraws %>% select("Population", "Asia Pacific", "Europe", "America"), level_plt, paste0(tail(levels(df_complete$EZ), 1), " by Geographic area"))

knitr::include_graphics(p_temp$p_temp_file)

Figures in S16

all_combination <- expand.grid(Experience = as.character(unique(df_complete$Experience)),
                               Country = as.character(unique(df_complete$Country))) %>% 
  as_tibble() %>% 
  left_join(propotion %>% select(Country, Continent), by = c("Country"))

mrp_pp <- posterior_linpred(mod_fit, transform = TRUE,
                            newdata = all_combination)

mrp_propotion <- all_combination %>% 
  left_join(propotion %>% 
              select(Country, Centers), 
            by = "Country")

mrp_pp <- mrp_pp %*% mrp_propotion$Centers / sum(mrp_propotion$Centers)

rawprop <- mean(df_complete$EZ == tail(levels(df_complete$EZ), 1))

p_temp <- plot_hist(mrp_pp[,1], ref = rawprop, paste0("MRP calibrated proportion of ", tail(levels(df_complete$EZ), 1)))

knitr::include_graphics(p_temp$p_temp_file)

SOZ defintion

df_complete <- df %>% 
  filter(SOZ != "Others",
         Country != "Unavailable",
         Experience != "Unavailable") %>% 
  mutate_if(is.factor, droplevels)

mod_fit <- brm(
  data = df_complete, family = bernoulli(link = "logit"),
  SOZ ~ 1 + (1 | Continent) + (1 | Country) + (1 | Experience),
  prior = c(prior(normal(0, 1.5), class = Intercept), 
            prior(normal(0, 1), class = sd)), 
  control = list(adapt_delta = 0.99),
  iter = 3500, warmup = 1000, chains = 4, cores = 4, seed = 20230918)

mod_fit

Model fitting results and diagnostics

 Family: bernoulli 
  Links: mu = logit 
Formula: SOZ ~ 1 + (1 | Continent) + (1 | Country) + (1 | Experience) 
   Data: df_complete (Number of observations: 296) 
  Draws: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
         total post-warmup draws = 10000

Group-Level Effects: 
~Continent (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.34      0.33     0.01     1.23 1.00     4694     4515

~Country (Number of levels: 34) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.54      0.26     0.07     1.10 1.00     2670     2792

~Experience (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.30      0.30     0.01     1.10 1.00     4955     5915

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     0.36      0.38    -0.42     1.12 1.00     5386     5498

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

fdraws <- as_draws_df(mod_fit) %>% 
  as_tibble() %>% 
  select(-starts_with(".")) %>% 
  select(-starts_with("lp")) %>% 
  summarise(Population = inv_logit_scaled(b_Intercept),
            `Low` = inv_logit_scaled(b_Intercept + `r_Experience[Low,Intercept]`),
            `Moderate` = inv_logit_scaled(b_Intercept + `r_Experience[Moderate,Intercept]`),
            `High` = inv_logit_scaled(b_Intercept + `r_Experience[High,Intercept]`)) 

level_plt <- c("Population", "Low", "Moderate", "High")

p_temp <- inference_plot_categorical(fdraws, level_plt, paste0(tail(levels(df_complete$SOZ), 1), " by Experience"))

knitr::include_graphics(p_temp$p_temp_file)

fdraws <- as_draws_df(mod_fit) %>% 
  as_tibble() %>% 
  select(-starts_with(".")) %>% 
  select(-starts_with("lp")) %>% 
  mutate(Population = inv_logit_scaled(b_Intercept),
         "Asia Pacific" = inv_logit_scaled(b_Intercept + `r_Continent[Asia.Pacific,Intercept]`),
         "Europe" = inv_logit_scaled(b_Intercept + `r_Continent[Europe,Intercept]`),
         "America" = inv_logit_scaled(b_Intercept + `r_Continent[America,Intercept]`)) 

fdraws <- fdraws %>% 
  select(-starts_with(c("b_", "r_", "sd_")))

level_plt <- c("Population", "Asia Pacific", "Europe", "America")

p_temp <- inference_plot_categorical(fdraws %>% select("Population", "Asia Pacific", "Europe", "America"), level_plt, paste0(tail(levels(df_complete$SOZ), 1), " by Geographic area"))

knitr::include_graphics(p_temp$p_temp_file)

Figures in S16

all_combination <- expand.grid(Experience = as.character(unique(df_complete$Experience)),
                               Country = as.character(unique(df_complete$Country))) %>% 
  as_tibble() %>% 
  left_join(propotion %>% select(Country, Continent), by = c("Country"))

mrp_pp <- posterior_linpred(mod_fit, transform = TRUE,
                            newdata = all_combination)

mrp_propotion <- all_combination %>% 
  left_join(propotion %>% 
              select(Country, Centers), 
            by = "Country")

mrp_pp <- mrp_pp %*% mrp_propotion$Centers / sum(mrp_propotion$Centers)

rawprop <- mean(df_complete$SOZ == tail(levels(df_complete$SOZ), 1))

p_temp <- plot_hist(mrp_pp[,1], ref = rawprop, paste0("MRP calibrated proportion of ", tail(levels(df_complete$SOZ), 1)))

knitr::include_graphics(p_temp$p_temp_file)

Other features

Figure 4a

df_tmp <- df %>% select(starts_with("feature"))
colnames(df_tmp) <- feature_names$name
df_tmp <- df_tmp %>% pivot_longer(everything(), values_to = "value", names_to = "feature") %>% 
                 mutate(feature = ifelse(feature == "Spike-ripple/-gamma", 
                                         "Spike-ripple/Spike-gamma activity", feature)) %>% 
                 group_by(feature) %>% summarize(prop_yes = mean(value == "Yes")) %>% 
                 arrange(prop_yes)

df_tmp$feature <- factor(df_tmp$feature, levels = c("Post-ictal spikes", "Post-ictal suppression",
                  "Spike-ripple/Spike-gamma activity", "Propagation patterns", "Ripples", 
                  "Spikes/slow waves",
                  "Pre-ictal spikes", "Near DC shift", "Fast ripples", 
                  "Electro-clinical correlations", "Stimulation", "Earliest seizure onset", 
                  "LVFA", "Seizure onset patterns"))

ggbarplot(df_tmp, "feature", "prop_yes", orientation = "horiz", fill = "grey90") + 
   scale_y_continuous(limits = c(0, 1), expand = c(0, 0.05), breaks = seq(0, 1, 0.1)) + 
   labs(y = "Response (%)")+
   theme(
        legend.position = "none",
        plot.margin = margin(5.5, 50, 5.5, 20, "pt"),
        axis.title.y = element_blank(),
        panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_blank(),
        panel.background = element_blank(),
        text = element_text(size = 38, face = "plain", color = "black"),
        axis.text = element_text(color = "black", size = 38),
        plot.title.position = "plot",
        plot.title = element_text(size = 36),
        plot.subtitle = element_text(size = 28)) -> p_temp
  
p_temp_file <- tempfile(tmpdir = "./fig", fileext = '.png')
agg_png(p_temp_file, width = 22, height = 12, units = "in", res = 300)
suppressWarnings(print(p_temp))
invisible(dev.off())

knitr::include_graphics(p_temp_file)

figure_list1 <- c()
figure_list2 <- c()
figure_list3 <- c()

df_complete <- df %>% 
  filter(Country != "Unavailable",
         Experience != "Unavailable") %>% 
  mutate_if(is.factor, droplevels)

all_combination <- expand.grid(Experience = as.character(unique(df_complete$Experience)),
                               Country = as.character(unique(df_complete$Country))) %>% 
  as_tibble() %>% 
  left_join(propotion %>% select(Country, Continent), by = c("Country"))

for (ifeature in seq(1, 14)){
  name_tmp <- feature_names[[ifeature, "name"]]
  idx_tmp <- feature_names[[ifeature, "index"]]
  
  mod_fit <- brm(
    data = df_complete %>% 
      select(Continent, Country, Experience, idx_tmp) %>% 
      rename(feature = idx_tmp), 
    family = bernoulli(link = "logit"),
    feature ~ 1 + (1 | Continent) + (1 | Country) + (1 | Experience),
    prior = c(prior(normal(0, 1.5), class = Intercept), 
              prior(normal(0, 1), class = sd)), 
    control = list(adapt_delta = 0.99),
    iter = 3500, warmup = 1000, chains = 4, cores = 4, seed = 20230918)
  
  # experience
  fdraws <- as_draws_df(mod_fit) %>% 
    as_tibble() %>% 
    select(-starts_with(".")) %>% 
    select(-starts_with("lp")) %>% 
    summarise(Population = inv_logit_scaled(b_Intercept),
            `Low` = inv_logit_scaled(b_Intercept + `r_Experience[Low,Intercept]`),
            `Moderate` = inv_logit_scaled(b_Intercept + `r_Experience[Moderate,Intercept]`),
            `High` = inv_logit_scaled(b_Intercept + `r_Experience[High,Intercept]`)) 
  
  level_plt <- c("Population", "Low", "Moderate", "High")
  
  p_temp <- inference_plot_categorical(fdraws, level_plt,
                                       paste0(" ", name_tmp, " by Experience"))
  

  figure_list1 <- c(figure_list1, p_temp$p_temp_file)
  
  # country
  fdraws <- as_draws_df(mod_fit) %>% 
    as_tibble() %>% 
    select(-starts_with(".")) %>% 
    select(-starts_with("lp")) %>% 
    mutate(Population = inv_logit_scaled(b_Intercept),
           "Asia Pacific" = inv_logit_scaled(b_Intercept + `r_Continent[Asia.Pacific,Intercept]`),
           "Europe" = inv_logit_scaled(b_Intercept + `r_Continent[Europe,Intercept]`),
           "America" = inv_logit_scaled(b_Intercept + `r_Continent[America,Intercept]`)) 
  
  for (country in propotion$Country) {
    if (propotion$Responses[propotion$Country == country] >= 4){
      fdraws[[country]] <- inv_logit_scaled(fdraws[["b_Intercept"]] + 
                                              fdraws[[paste0("r_Continent[", 
                                                             str_replace(propotion$Continent[propotion$Country == country], " ", "."), 
                                                             ",Intercept]")]] + 
                                              fdraws[[paste0("r_Country[", 
                                                             str_replace(country, " ", "."), 
                                                             ",Intercept]")]])
    }
    
  }  
  
  fdraws <- fdraws %>% 
    select(-starts_with(c("b_", "r_", "sd_")))
  
  level_plt <- c("Population", "Asia Pacific", "Europe", "America")
  
  p_temp <- inference_plot_categorical(fdraws %>% 
                                         select("Population", "Asia Pacific", "Europe", "America"), 
                                       level_plt,
                                       paste0(" ", name_tmp, 
                                              " by Geographic area"))
  
  figure_list2 <- c(figure_list2, p_temp$p_temp_file)
  
  # MRP
  mrp_pp <- posterior_linpred(mod_fit, transform = TRUE,
                              newdata = all_combination)
  
  mrp_propotion <- all_combination %>% 
    left_join(propotion %>% 
                select(Country, Centers), 
              by = "Country")
  
  mrp_pp <- mrp_pp %*% mrp_propotion$Centers / sum(mrp_propotion$Centers)
  
  rawprop <- mean(df_complete[[idx_tmp]] == "Yes")
  
  p_temp <- plot_hist(mrp_pp[,1], ref = rawprop, name_tmp)
  
  figure_list3 <- c(figure_list3, p_temp$p_temp_file)

}

Figures in S15

knitr::include_graphics(figure_list1[1])

knitr::include_graphics(figure_list2[1])

knitr::include_graphics(figure_list1[2])

knitr::include_graphics(figure_list2[2])

knitr::include_graphics(figure_list1[3])

knitr::include_graphics(figure_list2[3])

knitr::include_graphics(figure_list1[4])

knitr::include_graphics(figure_list2[4])

knitr::include_graphics(figure_list1[5])

knitr::include_graphics(figure_list2[5])

knitr::include_graphics(figure_list1[6])

knitr::include_graphics(figure_list2[6])

knitr::include_graphics(figure_list1[7])

knitr::include_graphics(figure_list2[7])

knitr::include_graphics(figure_list1[8])

knitr::include_graphics(figure_list2[8])

knitr::include_graphics(figure_list1[9])

knitr::include_graphics(figure_list2[9])

knitr::include_graphics(figure_list1[10])

knitr::include_graphics(figure_list2[10])

knitr::include_graphics(figure_list1[11])

knitr::include_graphics(figure_list2[11])

knitr::include_graphics(figure_list1[12])

knitr::include_graphics(figure_list2[12])

knitr::include_graphics(figure_list1[13])

knitr::include_graphics(figure_list2[13])

knitr::include_graphics(figure_list1[14])

knitr::include_graphics(figure_list2[14])

Figures in S16

knitr::include_graphics(figure_list3[1])

knitr::include_graphics(figure_list3[2])

knitr::include_graphics(figure_list3[3])

knitr::include_graphics(figure_list3[4])

knitr::include_graphics(figure_list3[5])

knitr::include_graphics(figure_list3[6])

knitr::include_graphics(figure_list3[7])

knitr::include_graphics(figure_list3[8])

knitr::include_graphics(figure_list3[9])

knitr::include_graphics(figure_list3[10])

knitr::include_graphics(figure_list3[11])

knitr::include_graphics(figure_list3[12])

knitr::include_graphics(figure_list3[13])

knitr::include_graphics(figure_list3[14])

License

The work is copyrighted by the ANPHY Lab and licensed under the BSD 3-clauselicense.

Copyright (c), 2024, John Thomas Email: j.thomas@duke.edu and Birgit Frauscher Email: birgit.frauscher@duke.edu. All rights reserved.

Code created by Zhengchen Cai, Email: zhengchen.cai@mcgill.ca

Computing environment

sessionInfo()

R version 4.1.2 (2021-11-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 22631)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252 
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    
system code page: 936

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggforce_0.3.3      moonBook_0.3.2     scales_1.2.1       ragg_1.2.1.9000   
 [5] webr_0.1.6         ggsci_2.9          ggpubr_0.4.0.999   ggstatsplot_0.12.1
 [9] brms_2.19.0        Rcpp_1.0.10        tidybayes_3.0.4    forcats_0.5.1     
[13] stringr_1.5.0      dplyr_1.1.3        purrr_1.0.2        readr_2.1.1       
[17] tidyr_1.3.0.9000   tibble_3.2.1       ggplot2_3.4.3      tidyverse_1.3.1   
[21] readxl_1.3.1      

loaded via a namespace (and not attached):
  [1] estimability_1.3        coda_0.19-4             bit64_4.0.5            
  [4] knitr_1.43              dygraphs_1.1.1.6        multcomp_1.4-18        
  [7] data.table_1.14.2       inline_0.3.19           generics_0.1.3         
 [10] cowplot_1.1.1           callr_3.7.3             TH.data_1.1-0          
 [13] correlation_0.8.4       bit_4.0.4               tzdb_0.2.0             
 [16] xml2_1.3.2              lubridate_1.8.0         httpuv_1.6.11          
 [19] StanHeaders_2.26.13     assertthat_0.2.1        xfun_0.39              
 [22] hms_1.1.1               ggdist_3.3.0            jquerylib_0.1.4        
 [25] bayesplot_1.8.1         evaluate_0.21           promises_1.2.0.1       
 [28] fansi_1.0.4             dbplyr_2.1.1            igraph_1.3.1           
 [31] DBI_1.1.2               tmvnsim_1.0-2           htmlwidgets_1.6.2      
 [34] tensorA_0.36.2          stats4_4.1.2            paletteer_1.4.0        
 [37] ellipsis_0.3.2          crosstalk_1.2.0         backports_1.4.1        
 [40] fontLiberation_0.1.0    V8_4.3.0                insight_0.19.10        
 [43] markdown_1.1            RcppParallel_5.1.4      fontBitstreamVera_0.1.1
 [46] vctrs_0.6.3             sjlabelled_1.1.8        abind_1.4-5            
 [49] cachem_1.0.8            withr_2.5.0             vroom_1.5.7            
 [52] vcd_1.4-11              checkmate_2.0.0         emmeans_1.7.1-1        
 [55] xts_0.13.1              prettyunits_1.1.1       mnormt_2.0.2           
 [58] ztable_0.2.3            crayon_1.5.2            crul_1.3               
 [61] labeling_0.4.2          pkgconfig_2.0.3         tweenr_1.0.2           
 [64] nlme_3.1-153            statsExpressions_1.5.2  rlang_1.1.1            
 [67] editData_0.1.8          lifecycle_1.0.3         miniUI_0.1.1.1         
 [70] MatrixModels_0.5-0      colourpicker_1.1.1      sandwich_3.0-1         
 [73] fontquiver_0.2.1        httpcode_0.3.0          modelr_0.1.8           
 [76] cellranger_1.1.0        distributional_0.3.2    polyclip_1.10-0        
 [79] lmtest_0.9-39           matrixStats_0.61.0      flextable_0.8.6        
 [82] datawizard_0.10.0       Matrix_1.4-0            loo_2.7.0.9000         
 [85] carData_3.0-4           boot_1.3-28             zoo_1.8-12             
 [88] reprex_2.0.1            base64enc_0.1-3         gamm4_0.2-6            
 [91] ggridges_0.5.3          processx_3.8.1          parameters_0.21.6      
 [94] shinystan_2.6.0         rstatix_0.7.0           ggsignif_0.6.4         
 [97] memoise_2.0.1           magrittr_2.0.3          plyr_1.8.7             
[100] threejs_0.3.3           compiler_4.1.2          rstantools_2.1.1       
[103] RColorBrewer_1.1-3      lme4_1.1-31             cli_3.6.1              
[106] pbapply_1.5-0           patchwork_1.1.1         ps_1.7.5               
[109] Brobdingnag_1.2-6       formatR_1.14            MASS_7.3-54            
[112] mgcv_1.8-38             tidyselect_1.2.0        stringi_1.7.12         
[115] textshaping_0.3.6       projpred_2.4.0.9000     yaml_2.3.7             
[118] askpass_1.1             svUnit_1.0.6            bridgesampling_1.1-2   
[121] grid_4.1.2              sass_0.4.6              tools_4.1.2            
[124] parallel_4.1.2          rio_0.5.29              rvg_0.3.2              
[127] rstudioapi_0.14.0-9000  uuid_1.0-3              foreign_0.8-81         
[130] gridExtra_2.3           devEMF_4.2              posterior_1.4.1        
[133] farver_2.1.0            digest_0.6.31           shiny_1.7.4            
[136] quadprog_1.5-8          gfonts_0.2.0            car_3.0-12             
[139] broom_0.7.12            BayesFactor_0.9.12-4.4  later_1.3.1            
[142] shinyWidgets_0.7.6      httr_1.4.2              gdtools_0.3.1          
[145] psych_2.1.9             effectsize_0.8.7        colorspace_2.0-2       
[148] rvest_1.0.2             fs_1.6.2                splines_4.1.2          
[151] rematch2_2.1.2          shinythemes_1.2.0       systemfonts_1.0.3      
[154] xtable_1.8-4            jsonlite_1.8.6          nloptr_1.2.2.3         
[157] rstan_2.26.13           zeallot_0.1.0           R6_2.5.1               
[160] pillar_1.9.0            htmltools_0.5.5         mime_0.12              
[163] glue_1.6.2              fastmap_1.1.1           minqa_1.2.4            
[166] DT_0.20                 codetools_0.2-18        pkgbuild_1.3.1         
[169] mvtnorm_1.1-3           utf8_1.2.3              lattice_0.20-45        
[172] bslib_0.5.0             arrayhelpers_1.1-0      curl_5.0.1             
[175] gtools_3.9.2            officer_0.6.0           zip_2.2.0              
[178] rrtable_0.3.0           shinyjs_2.0.0           openxlsx_4.2.5         
[181] openssl_2.0.6           survival_3.2-13         rmarkdown_2.22.2       
[184] munsell_0.5.0           sjmisc_2.8.9            haven_2.4.3            
[187] reshape2_1.4.4          gtable_0.3.1            bayestestR_0.13.2